{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "75c277ac",
   "metadata": {},
   "source": [
    "## Multivariate Linear Model - MultivariateLS\n",
    "\n",
    "This notebooks illustrates some features for the multivariate linear model estimated by least squares. \n",
    "The example is based on the UCLA stats example in https://stats.oarc.ucla.edu/stata/dae/multivariate-regression-analysis/ .\n",
    "\n",
    "The model assumes that a multivariate dependent variable depends linearly on the same set of explanatory variables.\n",
    "\n",
    "Y = X * B + u\n",
    "\n",
    "where  \n",
    "- the dependent variable (endog) `Y` has shape (nobs, k_endog),  \n",
    "- the matrix of explanatory variables including constant (exog) `X` has shape (nobs, k_exog), and\n",
    "- the parameter matrix `B` has shape (k_exog, k_endog), i.e. coefficients for explanatory variables in rows and dependent variables in columns.\n",
    "- the disturbance term `u` has the same shape as `Y`, (nobs, k_endog), and is assumed to have mean zero and to be uncorrelated with the exog `X`.\n",
    "\n",
    "Estimation is by least squares. The parameter estimates with common explanatory variables for each dependent variables corresponds to separate OLS estimates for each `endog`. The main advantage of the multivariate model is that we can make inference "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e644acb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "from statsmodels.base.model import LikelihoodModel\n",
    "from statsmodels.regression.linear_model import OLS\n",
    "from statsmodels.multivariate.manova import MANOVA\n",
    "from statsmodels.multivariate.multivariate_ols import MultivariateLS\n",
    "\n",
    "import statsmodels.multivariate.tests.results as path\n",
    "\n",
    "dir_path = os.path.dirname(os.path.abspath(path.__file__))\n",
    "csv_path = os.path.join(dir_path, \"mvreg.csv\")\n",
    "data_mvreg = pd.read_csv(csv_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e851e0ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_mvreg.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68ec2582",
   "metadata": {},
   "outputs": [],
   "source": [
    "formula = \"locus_of_control + self_concept + motivation ~ read + write + science + prog\"\n",
    "mod = MultivariateLS.from_formula(formula, data=data_mvreg)\n",
    "res = mod.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71687b33",
   "metadata": {},
   "source": [
    "### Multivariate hypothesis tests mv_test\n",
    "\n",
    "The `mv_test` method by default performs the hypothesis tests that each term in the formula has no effect on any of the dependent variables (`endog`). This is the same as the MANOVA test.  \n",
    "Note, MANOVA in statsmodels is implemented as test on coefficients in the multivariate model and is not restricted to categorical variables. In the current example, we have three continuous and one categorical explanatory variables, in addition to the constant.\n",
    "\n",
    "Consequently, using mv_test in MultivariateLS and in MANOVA produces the same results.\n",
    "However. MANOVA only provides the hypothesis tests as feature, while MultivariateLS provide the usual model results.\n",
    "\n",
    "More general versions of the mv_test are for hypothesis in the form `L B M = C`.\n",
    "Here `L` are restrictions corresponding to explanatory variables, `M` are restrictions corresponding to dependent (endog) variables and `C` is a matrix of constants for affine restrictions. See docstrings for details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cb94da2",
   "metadata": {},
   "outputs": [],
   "source": [
    "mvt = res.mv_test()\n",
    "mvt.summary_frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b999e805",
   "metadata": {},
   "outputs": [],
   "source": [
    "manova = MANOVA.from_formula(formula, data=data_mvreg)\n",
    "manova.mv_test().summary_frame"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff1664a9",
   "metadata": {},
   "source": [
    "The core multivariate regression results are displayed by the `summary` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "015cd62b",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fb028ed",
   "metadata": {},
   "source": [
    "The the standard results attributes for the parameter estimates like `params`, `bse`, `tvalues` and `pvalues`, are two dimensional arrays or dataframes with explanatory variables (`exog`) in rows and dependend (`endog`) variables in columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1295b73",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbb9042d",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.bse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9929394",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.pvalues"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd52087a",
   "metadata": {},
   "source": [
    "### General MV and Wald tests \n",
    "\n",
    "The multivariate linear model allows for multivariate test in the `L B M` form as well as standard wald tests on linear combination of parameters.  \n",
    "\n",
    "The multivariate tests are based on eigenvalues or trace of the matrices. Wald tests are standard test base on the flattened (stacked) parameter array and their covariance, hypothesis are of the form `R b = c` where `b` is the column stacked parameter array. The tests are asymptotically equivalent under the model assumptions but differ in small samples.\n",
    "\n",
    "The linear restriction can be defined either as hypothesis matrices (numpy arrays) or as strings or list of strings.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0a1d9f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "yn = res.model.endog_names\n",
    "xn = res.model.exog_names\n",
    "yn, xn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02982602",
   "metadata": {},
   "outputs": [],
   "source": [
    "# test for an individual coefficient\n",
    "\n",
    "mvt = res.mv_test(hypotheses=[(\"coef\", [\"science\"], [\"locus_of_control\"])])\n",
    "mvt.summary_frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b31d01db",
   "metadata": {},
   "outputs": [],
   "source": [
    "tt = res.t_test(\"ylocus_of_control_science\")\n",
    "tt, tt.pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c22a61e4",
   "metadata": {},
   "source": [
    "We can use either mv_test or wald_test for the joint hypothesis that an explanatory variable has no effect on any of the dependent variables, that is all coefficient for the explanatory variable are zero.\n",
    "\n",
    "In this example, the pvalues agree at 3 decimals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80b1e726",
   "metadata": {},
   "outputs": [],
   "source": [
    "wt = res.wald_test(\n",
    "    [\"ylocus_of_control_science\", \"yself_concept_science\", \"ymotivation_science\"],\n",
    "    scalar=True,\n",
    ")\n",
    "wt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12c0f058",
   "metadata": {},
   "outputs": [],
   "source": [
    "mvt = res.mv_test(hypotheses=[(\"science\", [\"science\"], yn)])\n",
    "mvt.summary_frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "335d0361",
   "metadata": {},
   "outputs": [],
   "source": [
    "# t_test provides a vectorized results for each of the simple hypotheses\n",
    "\n",
    "tt = res.t_test(\n",
    "    [\"ylocus_of_control_science\", \"yself_concept_science\", \"ymotivation_science\"]\n",
    ")\n",
    "tt, tt.pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "caf162f5",
   "metadata": {},
   "source": [
    "**Warning:** the naming pattern for the flattened parameter names used in `t_test` and `wald_test` might still change.\n",
    "\n",
    "The current pattern is `\"y{endog_name}_{exog_name}\"`.\n",
    "\n",
    "examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ca87db4",
   "metadata": {},
   "outputs": [],
   "source": [
    "[f\"y{endog_name}_{exog_name}\" for endog_name in yn for exog_name in [\"science\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6f7418d",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = [\n",
    "    f\"y{endog_name}_{exog_name}\"\n",
    "    for endog_name in yn\n",
    "    for exog_name in [\"prog[T.general]\", \"prog[T.vocational]\"]\n",
    "]\n",
    "c"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d0d1529",
   "metadata": {},
   "source": [
    "The previous restriction corresponds to the MANOVA type test that the categorical variable \"prog\" has no effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "960cb28e",
   "metadata": {},
   "outputs": [],
   "source": [
    "mant = manova.mv_test().summary_frame\n",
    "mant.loc[\"prog\"]  # [\"Pr > F\"].to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2b5ae70",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.wald_test(c, scalar=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab386919",
   "metadata": {},
   "source": [
    "**Note:** The degrees of freedom differ across hypothesis test methods.\n",
    "The model can be considered as a multivariate model with nobs=600 in this case, or as a stacked model with \n",
    "nobs_total = nobs * k_endog = 1800.\n",
    "\n",
    "\n",
    "For within endog restriction, inference is based on the same covariance of the parameter estimates in MultivariateLS and OLS. The degrees of freedom in a single output OLS are df_resid = 600 - 6 = 594. Using the same degrees of freedom in MultivariateLS preserves the equivalence for the analysis of each endog. Using the total df_resid for hypothesis tests would make them more liberal.\n",
    "\n",
    "Asymptotic inference based on normal and chisquare distribution (`use_t=False`) is not affected by how df_resid are defined.\n",
    "\n",
    "It is not yet decided whether there will be additional options to choose different degrees of freedom in the Wald tests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0854ceb",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.df_resid"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b5411b1",
   "metadata": {},
   "source": [
    "Both mv_test and wald_test can be used to test hypothesis on contrasts between coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8deb82c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = [\n",
    "    f\"y{endog_name}_prog[T.general] - y{endog_name}_prog[T.vocational]\"\n",
    "    for endog_name in yn\n",
    "]\n",
    "c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e0ee5d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.wald_test(c, scalar=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf983044",
   "metadata": {},
   "outputs": [],
   "source": [
    "mvt = res.mv_test(\n",
    "    hypotheses=[(\"diff_prog\", [\"prog[T.general] - prog[T.vocational]\"], yn)]\n",
    ")\n",
    "mvt.summary_frame"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1569e45",
   "metadata": {},
   "source": [
    "Example: hypothesis that coefficients are the same across endog equations.\n",
    "\n",
    "We can test that the difference between the parameters of the later two equation with the first equation are zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75ec5fe9",
   "metadata": {},
   "outputs": [],
   "source": [
    "mvt = res.mv_test(\n",
    "    hypotheses=[\n",
    "        (\n",
    "            \"diff_prog\",\n",
    "            xn,\n",
    "            [\"self_concept - locus_of_control\", \"motivation - locus_of_control\"],\n",
    "        )\n",
    "    ]\n",
    ")\n",
    "mvt.summary_frame"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a0b8193",
   "metadata": {},
   "source": [
    "In a similar hypothesis test, we can test that equation have the same slope coefficients but can have different constants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c8aa712",
   "metadata": {},
   "outputs": [],
   "source": [
    "xn[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04315196",
   "metadata": {},
   "outputs": [],
   "source": [
    "mvt = res.mv_test(\n",
    "    hypotheses=[\n",
    "        (\n",
    "            \"diff_prog\",\n",
    "            xn[1:],\n",
    "            [\"self_concept - locus_of_control\", \"motivation - locus_of_control\"],\n",
    "        )\n",
    "    ]\n",
    ")\n",
    "mvt.summary_frame"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c55852a",
   "metadata": {},
   "source": [
    "### Prediction\n",
    "\n",
    "\n",
    "The regression model and its results instance have methods for prediction and residuals.\n",
    "\n",
    "Note, because the parameter estimates are the same as in the OLS estimate for individual endog, the predictions will also be the same between the MultivariateLS model and the set of individual OLS models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba92e621",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.resid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa8885f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.predict()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bda82cf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.predict(data_mvreg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b0fec2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.fittedvalues"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa6d30de",
   "metadata": {},
   "source": [
    "The predict methods can take user provided data for the explanatory variables, but currently are not able to automatically create sets of explanatory variables for interesting effects.\n",
    "\n",
    "In the following, we construct at dataframe that we can use to predict the conditional expectation of the dependent variables for each level of the categorical variable \"prog\" at the sample means of the continuous variables. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ae189fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_exog = data_mvreg[[\"read\", \"write\", \"science\", \"prog\"]]\n",
    "\n",
    "ex = pd.DataFrame(data_exog[\"prog\"].unique(), columns=[\"prog\"])\n",
    "mean_ex = data_mvreg[[\"read\", \"write\", \"science\"]].mean()\n",
    "\n",
    "ex.loc[:, [\"read\", \"write\", \"science\"]] = mean_ex.values\n",
    "ex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "773618e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "pred = res.predict(ex)\n",
    "\n",
    "pred.index = ex[\"prog\"]\n",
    "pred.columns = res.fittedvalues.columns\n",
    "print(\"predicted mean by 'prog':\")\n",
    "pred"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69a05b73",
   "metadata": {},
   "source": [
    "## Outlier-Influence\n",
    "\n",
    "This is currently in draft version.  \n",
    "`resid_distance` is a one dimensional residual measure based on Mahalanobis distance for each sample observation. \n",
    "The hat matrix in the MultivariateLS model is the same as in OLS, the diagonal of the hat matrix is temporarily attached as `results._hat_matrix_diag`.\n",
    "\n",
    "Note, individual components of the multivariate dependent variable can be analyzed with OLS and are available in the corresponding post-estimation methods like `OLSInfluence`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d69a12c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.resid_distance[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f885722",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.cov_resid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "542b5eb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d21fc8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(res.resid_distance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adff83a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(res._hat_matrix_diag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b56d845f",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(res._hat_matrix_diag, res.resid_distance, \".\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
